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Temporal and spatial dynamics of scaling-specific 
features of a gene regulatory network in Drosophila 


Honggang Wu!23, Manu‘, Renjie Jiao”? & Jun Ma6 


A widely appreciated aspect of developmental robustness is pattern formation in proportion 
to size. But how such scaling features emerge dynamically remains poorly understood. Here 
we generate a data set of the expression profiles of six gap genes in Drosophila melanogaster 
embryos that differ significantly in size. Expression patterns exhibit size-dependent dynamics 
both spatially and temporally. We uncover a dynamic emergence of under-scaling in the 
posterior, accompanied by reduced expression levels of gap genes near the middle of large 
embryos. Simulation results show that a size-dependent Bicoid gradient input can lead to 
reduced Krtippel expression that can have long-range and dynamic effects on gap gene 
expression in the posterior. Thus, for emergence of scaled patterns, the entire embryo may be 
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viewed as a single unified dynamic system where maternally derived size-dependent infor- 
mation interpreted locally can be propagated in space and time as governed by the dynamics 
of a gene regulatory network. 
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individual’s body size is an intriguing feature of animal 

development!~*. There are two aspects directly relevant to 
this developmental scaling problem—scaled tissue specification 
and coordinated tissue expansion®®. In many developmental 
systems, the events controlling these two aspects are connected 
both temporally and mechanistically’~!°, and recent studies have 
investigated mechanisms controlling the formation and action of 
morphogen gradients whose profiles are scaled with the 
expanding size of a tissue!!~4, In other developmental systems, 
tissue specification and tissue expansion can take place in a 
temporally sequential manner®. For example, formation of scaled 
patterns in the chick and mouse neural tube is controlled in two 
sequential phases, morphogen-induced progenitor specification 
followed by cell-type-specific growth!*. In Drosophila embryos, 
patterning along the anterior-posterior (AP) axis can be viewed 
purely as a patterning problem because the physical size of the 
patterning system, the embryo, is predetermined at an earlier 
stage of the life cycle (that is, oogenesis) and does not change 
during the time of pattern formation®!®. This system thus 
provides a unique window to probing how scaled patterns emerge 
from the dynamic operation of a gene regulatory network without 
the entanglement of growth. 

Gap genes are situated at the top of the zygotic regulatory 
hierarchy for instructing the AP patterning outcome in 
Drosophila!’-*+, These genes respond initially to maternal 
inputs such as the Bicoid (Bcd) morphogen gradient and, 
subsequently, rely on extensive cross-regulatory mechanisms 
that form an intricately wired regulatory network!?°-?7, 
Previous studies have generated useful insights into the 
regulatory mechanisms controlling the spatial and temporal 
dynamics of gap gene expression?™?®-30, But how these dynamics 
are related to, and impacted by, embryo size remains poorly 
understood. In addition, while the relationship between the 
scaling properties of the Bcd gradient input and the expression of 
the gap gene hunchback (hb) has been studied quantitatively®>', 
whether (and if so how) this maternal input may impact the 
dynamic operation of the entire gap gene network in relation to 
embryo size remains to be resolved. 

The concept of gene regulatory networks provides a mechan- 
istic view of how the precise and reproducible patterning and cell 
fate specification are controlled”*. In addition to the scaling 
problem in embryonic pattern formation that we investigate in 
the current study, developmental robustness can also be 
understood in the context of other types of variability that a 
developmental system must face, such as molecular noise*2~*, 
environmental fluctuations!®7°3” and genetic variations*®?. 
While the term robustness (also often referred to as 
canalization) projects an effective image of developmental 
systems producing reliable outcomes in spite of the various 
uncertainties, robustness may not be viewed as an absolute 
concept. There is a ‘sensitive’ side of developmental systems that 
is exemplified by concentration-dependent actions of regulatory 
proteins?+!32,49.41 or human diseases caused by either having an 
extra copy or mutating a single copy of a regulatory gene*”. 
Thus experimental and theoretical studies that probe molecular 
origins and limits of robustness in well-defined systems can lead 
to fundamental insights into how developmental programs are 
controlled. With respect to the scaling problem in Drosophila AP 
patterning that we investigate here, accumulating evidence 
suggests that scaling is imperfect®?®39444, but the underlying 
molecular mechanisms are not fully resolved, particularly for the 
posterior part of the embryo®. Since gap genes are the earliest 
responders to maternal inputs and form a well-characterized 
regulatory network, systematically evaluating their expression 
properties in relation to embryo size will likely lead to new 
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insights into how size-dependent (that is, scaling-specific) 
features of a patterning system emerge. Such new insights 
can strengthen our knowledge of the concept of develop- 
mental robustness in terms of its underlying molecular origins 
and limits. 

Here we perform quantitative studies in Drosophila embryos 
that differ significantly in egg length and permit reliable 
extraction of scaling-specific information. We quantify the 
expression properties of six gap genes, orthodenticle (otd), hb, 
giant (gt), Kriippel (Kr), knirps (kni) and tailless (tll) in these 
embryos. Our results reveal that scaling-specific features of gap 
gene expression are highly dynamic in both space and time. We 
show that gap expression boundaries in the anterior exhibit a 
general over-scaling (over-scaling is defined as an excessive 
compensation for embryo length to give rise to a more posteriorly 
located boundary in large embryos relative to small embryos), 
which is consistent with the scaling properties of the Bcd gradient 
input®?!46, The expression boundaries in the posterior trunk 
region of the embryo undergo a spatially and temporally 
concerted dynamic under-scaling (under-scaling is defined as 
an insufficient compensation for embryo length to give rise to a 
more anteriorly located boundary in large embryos relative to 
small embryos). We perform simulation studies to probe the 
origins and impacts of an experimentally detected divergence in 
gap gene expression levels between large and small embryos. Our 
results support a view that the divergence in Kr expression level is 
a ‘localized’ scaling-specific feature that can originate from the 
size-dependent properties of the Bcd gradient and exert long- 
range and dynamic effects on gap gene expression behaviour in 
the posterior. Our study suggests that the molecular interactions 
(of the gap gene network) that control gene expression behaviour 
within individual cells (nuclei) also mediate the propagation of 
size-dependent ‘information’ in space and time throughout the 
embryo. Thus, the entire embryo may be viewed as a single 
unified dynamic system, in which a ‘localized’ interpretation of a 
size-dependent morphogen input can impact the dynamics and 
outcome of the patterning landscape as a whole. Such a ‘holistic’ 
view may be of general importance to studies of problems in 
developmental biology. 


Results 
Experimental design and data set. To facilitate experimental 
investigation of scaling-specific properties of AP patterning, 
we took advantage of two Drosophila inbred lines that had 
been artificially selected to lay large and small eggs*46 
(hereafter referred to as the large-egg and small-egg line, 
respectively). Embryos from these two lines are referred to as 
the large and small embryos, respectively. Their mean lengths, 
482.0 +20.9um (mean s.d.) and 408.6 +16.8um, differ by 
~15% (P=2.94x 1076; length measurements are from T3 
embryos for hb FISH (fluorescence in situ hybridization); see 
Supplementary Table 1 for n and below for further details). We 
used whole mount mRNA FISH to analyse the expression 
patterns of six gap genes that form boundaries distributed across 
the entire length of the embryo. Our experimental and imaging 
procedures were designed to avoid any nonlinear steps of 
signal amplification or adjustments. Under our experimental 
and analytical framework*!, the extracted intensities of the 
fluorescent signals preserve a linear relationship with the 
cytoplasmic mRNA concentrations (in arbitrary units, a.u.). To 
permit direct comparisons, both sets of embryos were treated 
side-by-side at both experimental and imaging steps. 

To facilitate the reconstruction of the temporal evolution of 
gene expression patterns, we sorted embryos into well-resolved 
temporal classes using morphological markers. This study focuses 
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primarily on nuclear cycles (nc) 13 and 14, during which gap gene 
transcripts have begun to accumulate to levels suitable for reliable 
quantitative measurements. This is also the time during which 
gap genes have begun to cross-regulate among themselves in an 
intensifying way!’. For embryos at ncl4, we divided them into 
nine time classes in a manner to ensure that each time class had a 
sufficient sample size for each and every gene analysed 
(Supplementary Table 1; see Methods for estimated time 
durations). In total, the data set generated in this study consists 
of data from 2,202 embryos, with 20 expression boundaries 
spanning across the entire AP length and across 10 time internals. 
Our data set permitted us to extract quantitative information 
about gap gene expression with regards to (1) their spatial 
features (that is, expression domains and boundary positions at a 
given time), (2) their temporal dynamics (that is, movements of 
domain boundaries and threshold-crossing positions as a 
function of time) and (3) their scaling-specific features (that is, 
relationships between the above mentioned properties and 
embryo size). 


Spatial registry and movements of gap expression boundaries. 
Figure 1a,b shows superimposed profiles of gap gene expression 
as a function of relative AP position in large and small 
embryos, respectively. Throughout this report, relative positions 
along the AP axis are expressed as values of fractional embryo 
length, ¢=x/L, with 0 denoting the anterior pole and 1 
the posterior pole. Our results shown in Fig. la,b yield the 
following spatial registry of expression boundaries for both large 
and small embryos along the AP axis: knil>otd1>tlll> 


hb1/gt3 > tll2 > otd2 > gt4 > hb2/Kr1 > Kr2/kni2 > kni3/gt5 > gt6/ 
hb3>tll3>hb4 (see Supplementary Fig. 1 for boundary 
nomenclature). Two additional boundaries, gtl and gt2, which 
emerge at later times in ncl4, can be included in this spatial 
registry that would remain identical between the two sets of 
embryos (Supplementary Fig. 1). These results show that, despite 
the considerable difference between their average lengths, 
embryos from the two inbred lines have an identical spatial 
registry of the expression boundaries of gap genes. They indicate 
that, at this initial level of analysis and consistent with the well- 
acknowledged framework of developmental robustness, expres- 
sion patterns of gap genes are generally robust to variations in 
embryo lengths. 

To analyse quantitatively the expression boundaries of 
gap genes and characterize their temporal movements 
during development, we measured their relative positions in 
individual embryos (Fig. 1; Methods). Figure 1c,d shows the mean 
positions of each of the analysed boundaries in the large and 
small embryos, respectively, at each of the time classes 
(see Supplementary Fig. 2 for s.d.). These results provide a visual 
representation summarizing how the gap expression boundaries 
move along the AP axis as a function of time. We can see readily 
that the boundaries located in the anterior and posterior parts of 
the embryos tend to exhibit movements toward mid-embryo (for 
example, ofd2 and hb3). This is in contrast to boundaries near 
mid-embryo (for example, gt4, hb2 and Kr1) or near the anterior 
cap (knil), whose positions remain relatively stable as a function 
of time. The overall dynamic movements of gap expression 
boundaries in embryos from the two inbred lines are generally 
similar to those in other characterized wild-type strains**-°°. 
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Figure 1 | Spatial registry and positional dynamics of gap gene expression boundaries. (a,b) Shown are expression profiles of the indicated six gap genes 
in large (a) or small (b) embryos at the peak time (in expression) of ncl4 (T2 and T3 combined) as a function of relative AP position €=x/L (O and 1 
denote anterior and posterior, respectively). For each gene, the intensity data extracted from an embryo along the AP axis were background-subtracted 
without any further adjustments. The mean and s.d. (shown as error bars) of the intensity data of the indicated genes are plotted for each of the two groups 
of embryos. (c,d) Shown are mean positions of gap expression boundaries in large (c) or small (d) embryos at the indicated time classes. Anterior and 
posterior boundaries are shown as left- and right-pointing triangles, respectively. Note that not all boundaries are present at all times (see Methods). 
Colour code from a to d: red, hb; magenta, gt; green, Kr; cyan, kni; blue, otd; black, tll. See Supplementary Table 1 for n. 
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Location-dependent scaling of gap expression boundaries. To 
evaluate AP patterning properties in relation to embryo length, 
we compared directly the positions of gap expression boundaries 
in large and small embryos. Our analysis reveals measurable 
differences between large and small embryos in either the 
boundary positions at a given time or their dynamic movements 
over time (Supplementary Fig. 2, and see legend for discussions of 
specific examples). In addition, the moving span of a boundary 
(that is, the shift of a boundary’s relative AP position during a 
period of time) can also differ between large and small embryos 
(Supplementary Fig. 2 legend). 

To quantify the relationship between AP patterning properties 
and embryo size, we employed a parameter, the scaling coefficient 
S, which was measured by pooling large and small embryos for a 
given time class (Methods; Supplementary Fig. 3). Under our 
current definition’, S=0 denotes perfect scaling, that is, a 
boundary position (expressed as a relative AP position č) exhibits 
a sample-wide insensitivity to variations in embryo length“. 
When S<0 or S>0, the boundary is under- or over-scaled, 
respectively, indicating that an increase in embryo length is either 
inadequately (under) or excessively (overly) compensated at the 
boundary position of a gene expression domain. Therefore, S 
quantifies exclusively boundary position variations derived from 
imperfect scaling despite the existence of other inevitable sources 
of noise that can also contribute to boundary variations (see 
Supplementary Fig. 4 for results supporting the conclusion that 
scaling contributes to the overall precision of gap expression 
boundaries). 

Figure 2a-t shows the S$ values for individual expression 
boundaries plotted as a function of time, with the 95% confidence 
intervals (CI) shown as shaded bands. These profiles reveal 
diverging scaling properties in the anterior and posterior parts 
of the embryo. In particular, while the gene expression 
boundaries in the anterior have a tendency of being over-scaled 
(that is, S>0), those in the posterior tend to be under-scaled (that 
is, S<0). In addition, expression boundaries located near the 
middle of the embryo (for example, gt4) or close to either pole 
(for example, gtl, hb4) have a tendency of exhibiting a nearly 
perfect scaling, particularly at later times of ncl4. To provide a 
visual representation of these scaling trends, we generated an 
illustrative, superimposed plot of the S profiles for gap expression 
boundaries (Fig. 2u). Here the over- and under-scaled boundaries 
become readily perceivable. These results show that gap 
expression boundaries located in different parts of the embryo 
have different scaling-specific properties, that is, scaling is not a 
spatially uniform feature of the gap gene network. 


Spatially and temporally concerted posterior under-scaling. A 
noted feature of the results shown in Fig. 2 is that, for expression 
boundaries in the posterior trunk region of the embryo (for 
example, kni2, kni3 and gt6), the degree of under-scaling is the 
greatest around mid-ncl4 (that is, for kni3 at T5, 
S= — 0.19 + 0.06; error bar is 95% CI). At earlier times, many of 
these boundaries are not under-scaled, with some even exhibiting 
over-scaling. An example is gt6, which has the highest degree 
of over-scaling at T1 (S=0.33 +0.18). But this over-scaling 
vanishes precipitously (S= 0.02 + 0.13 at T2) to become under- 
scaled as embryos progress further into ncl4 (for example, 
S= — 0.14 Ł 0.07 at T4). Figure 3 shows S profiles for expression 
boundaries detected along the AP axis at individual time classes 
(see Supplementary Fig. 5 for a three-dimensional illustration). 
Here, under-scaling of gap expression boundaries in the posterior 
trunk region of the embryo at mid-to-late times of nc14 exhibits 
itself as one of the most visible features. This feature can be 
further appreciated when evaluated against the S profiles at earlier 


times (nc13, T1 and T2) before the emergence of under-scaling 
for these boundaries. 

Another feature displayed in Fig. 3 is a ‘wave-like’ under- 
scaling that emerges near mid-embryo and propagates toward the 
posterior. Specifically, a sudden dip (S amplitude trough) 
develops around middle embryo at T2. This S amplitude trough 
moves towards the posterior to cause a broader swath of the 
embryo become under-scaled. By the time of T6, the S amplitude 
trough has reached a location near č = ~ 0.8, while at the same 
time the region near mid-embryo has been fully recovered. The S 
amplitude trough never reaches the posterior pole region marked 
by the hb4 boundary (see also Fig. 2t), and portions of the 
posterior trunk region never recover fully. Together, our analyses 
of the S profiles show that the scaling-specific features of the AP 
patterning network emerge dynamically with regard to both space 
and time. 


Large embryos have reduced levels of gap expression in the 
middle. The location- and time-dependent scaling properties of 
gap expression boundaries described above represent new 
discoveries that have not been documented previously. They 
suggest that, despite extensive studies, our current knowledge 
about how this regulatory network operates remains incomplete. 
To probe potential mechanisms, we gauged the effects of embryo 
size on spatiotemporal patterning by analysing the alterations to 
gap expression domain widths in time (see Fig. 4a for a cartoon 
depicting a moving expression domain that has a ‘leading’ 
boundary and a ‘trailing’ boundary in a shifting landscape). Here 
we measured the widths of each expression domain (expressed as 
values relative to embryo length) in the large and small embryos 
to obtain a difference, AW (Methods; Supplementary Fig. 6). 
Figure 4b-d shows that AW for Kr, kni and gt domains in the 
posterior trunk region declines as a function of time, suggesting 
that these domains experience a greater degree of contraction in 
large embryos relative to small embryos. This effect for the Kr 
domain is most striking, where an initial positive AW becomes a 
negative value at later times of ncl4 (Fig. 4b). Figure 4e,f shows 
that, for expression domains closer to the posterior pole (that is, 
the posterior hb and tll domains), AW remains generally a 
positive value throughout ncl4 (Supplementary Fig. 6). 

To further probe potential mechanisms or events associated 
with the emergence of under-scaling in the posterior trunk 
region, we compared the expression levels of gap domains in large 
and small embryos (see Methods). Figure 5a shows that, relative 
to small embryos, large embryos exhibit properties of either 
decreased or progressively decreasing expression levels in the 
posterior trunk region. In particular, throughout nc14, Kr level is 
persistently lower in large embryos than in small embryos 
(Fig. 5a; Supplementary Fig. 1). While the expression levels of 
both the kni and hb domains at earlier times of ncl4 are 
comparable between large and small embryos, their levels in large 
embryos progressively become lower relative to small embryos 
(Fig. 5a). The posterior gf domain in large embryos exhibits a 
particularly dynamic decrease in expression output, with a higher 
level at earlier times but a lower level at later times (Fig. 5a). The 
fact that large embryos do not have universally (gene- and time- 
wise) low expression levels also rules out systematic experimental 
errors as the origin of differential expression levels. Together, the 
dynamic differences in gap domain widths (Fig. 4) and expression 
levels (Fig. 5a) suggest that, relative to small embryos, large 
embryos tend to have a weaker gene expression output 
particularly near the middle at the Kr domain. 


Evaluating potential originators of posterior under-scaling. To 
further evaluate which part of the embryo may play an originator 
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Figure 2 | Temporal evolution of scaling coefficient of gap expression boundaries. (a-t) Shown are scaling coefficient S profiles for individual gap 

expression boundaries at the indicated times. Panels are ranked according to the AP registry of the boundaries (top left most anterior, bottom right most 
posterior). A dotted line denotes S= 0 and shaded bands are 95% CI of the linear regression slope (Supplementary Fig. 3). (u) Superimposed plot of S 
profiles for all the boundaries from a to t (except the fast-decaying hb1), with two colours denoting boundaries located in the anterior (dark red) or the 


posterior trunk region (purple). 


role in the emergence of posterior under-scaling, we compared 
the positions at which a mean expression profile crosses given 
thresholds in large and small embryos. Unlike measurements of 
boundary position or domain width, which utilize normalized 
intensity data (Fig. 4a; Methods), threshold-crossing analysis was 
based on unnormalized data, thus also taking into account the 
expression level of a gene at a time. Figure 5b shows scatter plots 
of threshold-crossing positions of the trailing boundaries in the 
posterior trunk region as well as the leading #113 boundary 
(Supplementary Fig. 7). Here we focus on time classes T2-T4, 
during which the observed under-scaling begins to emerge (Figs 2 
and 3; Supplementary Fig. 5). 

Our analysis of threshold-crossing positions of gap expression 
profiles leads to the following findings. First, when differences 
between large and small embryos are considered, the trailing Kr 
boundary (that is, the posterior boundary of the central Kr 


domain) moves in a manner that is earliest in time and largest in 
moving span (see above for definition) among all the trailing 
boundaries in the posterior trunk region (Fig. 5b). Second, the 
threshold-crossing positional differences (between large and small 
embryos) for leading boundaries are generally less than trailing 
boundaries (Supplementary Fig. 7). Third and importantly, the 
threshold-crossing positions of the tll3 boundary do not expand 
appreciably more from the posterior pole in large embryos than 
in small embryos (Fig. 5b). These results argue against the 
posterior pole as the originator of under-scaling. 


Long-range and dynamic impacts of a lowered Kr level. Our 
results described thus far suggest that the size-dependent 
divergence in gap expression levels near mid-embryo (Fig. 5a) is 
itself a scaling-specific feature that could potentially account for 
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Figure 3 | Scaling coefficient profiles as a function of AP position. Shown are measured S profiles of gap expression boundaries plotted as a function of 
AP position at each of the 10 time classes, with blue bands representing 95% CI of S. The two dotted lines denote mid-embryo (relative AP position 


€=0.5) and perfect scaling (S= 0), respectively. 


some of the observed patterning behaviour in the posterior trunk 
region. To directly evaluate this possibility, we performed simu- 
lation studies based on a gene circuit model of the gap gene 
network that correctly simulates dynamic domain boundary 
shifts as well as the effects of embryo-to-embryo Bcd 
variation!?*”, Here we specifically reduced the maximal 
synthesis rates of gap proteins without any other adjustments 
to the model!4 (Methods). In essence, our analysis was specifically 
designed to simulate how a reduction in gap protein levels—as a 
given—affects the behaviour of the operation of the gap gene 
network. In our analysis, we reduced the maximal synthesis rates 
of one protein at a time (Kr, Kni, Gt or Hb) or reduced all four 
simultaneously. 

Our results are shown in Supplementary Fig. 8a-e, where 
protein expression profiles of different gap genes are colour-coded 
and plotted as a function of relative AP position. They reveal the 
following findings. First, individually reducing the synthesis rate of 
Hb or Gt in the embryo only had modest and locally restricted 
effects on gap gene expression patterns in the posterior 
(Supplementary Fig. 8a,b). Second, simultaneously reducing the 
synthesis rates of all four proteins did not cause a strong alteration 
in expression boundaries (Supplementary Fig. 8c). Third and 


importantly, individually reducing synthesis rate of either Kr or 
Kni led to long-range effects on gene expression boundaries in the 
posterior (Supplementary Fig. 8d,g). However, unlike Kr level 
reduction that could recapitulate qualitatively the anterior shifts of 
several expression boundaries (for example, Gt5, Gt6 and Hb3), 
reducing Kni level caused them to shift posteriorly (compare 
panels e and h in Supplementary Fig. 8). 

Reducing Kr level on its own can dynamically affect the 
expression levels of other gap proteins in the posterior. While the 
Kni expression level is lowered throughout nc14, Gt expression is 
higher at earlier times only to become lower at later times 
(Supplementary Fig. 8f). This simulated dynamic property is 
qualitatively very similar to the experimentally observed mRNA 
expression behaviour of these two genes in the posterior trunk 
region (Fig. 5a). Furthermore, reducing the synthesis rate of Kr 
can lower the Kni level but not the other way around 
(Supplementary Fig. 8fi). Together, these results suggest that 
the experimentally observed Kr mRNA level divergence between 
large and small embryos is a scaling-specific feature that on its 
own can exert long-range and dynamic effects on gap gene 
expression behaviour in the posterior (Supplementary Fig. 8 
legend). 
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Figure 4 | Dynamic changes in gap expression domain widths. (a) Shown is a schematic diagram of an expression domain illustrating the definition of its 
width (W), the direction of its movement during nc14, and a depiction of its leading and trailing boundaries. (b-f) Shown are differences in domain widths 
(expressed as values relative to embryo length) between large and small embryos as a function of time. Shaded band for each panel represents s.d. of AW. 
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Figure 5 | Large embryos exhibit weaker or weakening gap expression levels in the posterior. (a) Gap gene expression levels in the posterior trunk 
region at T2-T5. For each embryo, the expression level is estimated as the mean across the three highest intensity bins of a domain. Blue and red denote 
large and small embryos, respectively. Error bars are s.d. Results of Student's t-tests are given, with n.s., * and ** denoting not significant, P< 0.05 and 
P< 0.01, respectively. (b) Shown are AP positions at which a mean intensity profile crosses given thresholds. The trailing boundaries for Kr, kni, gt and hb as 
well as the leading boundary of posterior tll domain are shown. Measured boundary positions of the large and small embryos are indicated at the bottom of 
each panel for reference. Blue, large embryos; red, small embryos. See Supplementary Table 1 for n. 
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Bcd input properties can impact gap gene network through Kr. 
The long-range effects of Kr protein level reduction documented 
by our simulation studies raised an important mechanistic 
question. If the experimentally observed divergence in Kr mRNA 
level (Fig. 5a) was indeed responsible, at least partially, for some 
of the scaling-specific features in the posterior trunk region, what 
properties might have led to such a divergence? To address this 
question, we turned our attention to Bcd (refs 6,31,46). One 
important property of the Bcd gradient is that the large-embryo 
profile intersects with the small-embryo profile at a critical 
position®, such that large embryos have higher concentrations 
anterior to the critical position and lower concentrations 
posterior to it>*°. The critical position lies anterior to the Kr 
domain, so that large embryos have lower Bcd concentrations at 
the location of Kr expression than small embryos*®. Bcd is known 
to activate Kr reporter constructs’, a result supported by in 
silico reconstructions of the gap gene network from expression 
data!™2550, A positive role for Bcd in Kr transcription suggests 
that Bcd gradient properties may contribute to (or be responsible 
for) the experimentally observed Kr mRNA level divergence 
between large and small embryos. To directly test this possibility, 
we performed simulation studies where we adjusted the 
amplitude and the relative length scale of the Bcd gradient in 
ways that are consistent with experimentally observed gradient 
properties in large embryos. In essence, we compared two 
otherwise identical embryo systems except for a difference in the 
Bcd gradient input that is designed to mimic a difference in 
embryo size. 

Our simulation results (Fig. 6; Supplementary Fig. 9) show that 
the scaling properties of the Bcd input can cause Kr level 
reduction in the simulated large embryos. Importantly, as we 
moved the critical position further anterior and away from the Kr 
expression domain, Kr level became progressively lower in the 
simulated large embryos (Supplementary Fig. 9). This documents 
a causal effect. In addition, the impact of the scaling properties of 
the Bcd gradient is not restricted to Kr expression level in the 
simulated large embryos. Other gap proteins in the posterior 
trunk region also exhibit reduced levels (Fig. 6j; Supplementary 
Fig. 9c,fi). Importantly, while Kr has a lower level in the 
simulated large embryos throughout ncl4, other gap proteins 
(most notably Kni and Gt) exhibit dynamic weakening in a 
manner that resembles the experimentally observed behaviour 
(Fig. 5a). 

The expression boundaries of gap genes in the posterior trunk 
region of the simulated large embryos also exhibit differences 
relative to their small counterparts in their spatiotemporal 
dynamics. The properties of two boundaries, Kni2 and Gt5, are 
particularly noteworthy. These two boundaries in large embryos 
exhibit an anterior shift relative to small embryos (that is, they are 
under-scaled) at early times, but such an under-scaling becomes 
‘recovered’ at later times. These simulated dynamic properties are 
reminiscent of the ‘wave-like’ initiation and recovery of under- 
scaling observed experimentally (Fig. 3; Supplementary Fig. 4). 
Together these simulation studies provide important mechanistic 
insights into how the regulatory logic of the gap gene network 
can permit a ‘localized’ scaling-specific feature, namely Kr 
mRNA level divergence near mid-embryo, to exert its long-range 
effects dynamically. In addition, under the framework of the gene 
circuit model!, this ‘localized’ scaling-specific feature near mid- 
embryo can originate from the scaling properties of the Bcd 
gradient. 


Discussion 
Anterior-posterior patterning in Drosophila embryos represents 
an excellent experimental paradigm for studying regulatory 


network operation and developmental robustness. Despite 
extensive studies, it is poorly understood how the AP patterning 
network adapts its operation to the size of an embryo in achieving 
developmental robustness. Here we perform experimental and 
simulation studies to gain mechanistic insights into this question. 
We use two inbred lines that together offer an enhanced spread in 
embryo length within a species to allow an efficient extraction of 
scaling-specific features of gap expression properties. We show 
that scaling is neither absolute in space nor constant in time. We 
uncover a spatially and temporally concerted dynamic emergence 
of under-scaling of gap expression boundaries in the posterior 
trunk region, in general agreement with even-skipped (eve) 
expression behaviour in this part of the embryo°. An important 
design feature of our experimental data set, which is the first of its 
kind, is to permit direct comparisons of gap expression levels 
in relation to embryo size. Our experimental and simulation 
results support a view that the dynamic emergence of 
under-scaling in the posterior trunk originates ‘locally’ as a 
result of lowered expression levels of gap genes near the middle of 
large embryos. 

Our results show that gap expression boundaries located in the 
anterior are generally over-scaled (Fig. 2), which is consistent 
with the scaling properties of the Bcd gradient input in this part 
of the embryo™®. In the anterior, the tll2 and otd2 boundaries 
are also sensitive to input from the terminal system?!*?, 
These boundaries exhibit posterior movements (Fig. 1c,d; 
Supplementary Fig. 2) that coincide with an intensification of 
over-scaling as embryos progress into nc14 (Fig. 2). These results 
are consistent with the suggestion that Bcd and terminal inputs 
interact to specify patterning properties in the anterior of the 
embryo*!-?, But the scaling properties of the Bcd gradient are 
inadequate to explain the scaling-specific features in the 
posterior®. It thus remained an open question as to how such 
features might emerge dynamically and how they might be 
controlled. The extensive interactions between gap genes in the 
posterior trunk region!™?526:53-55 suggest that regulatory 
networks and/or maternal inputs other than Bcd may play 
critical roles in patterning this part of the embryo. But our 
observed dynamic emergence of under-scaling is a new finding 
and, thus, a fresh conceptual framework is required to gain 
mechanistic insights. Our simulation results suggest that the 
scaling properties of the Bcd gradient input can lead to the 
experimentally observed Kr mRNA level divergence between large 
and small embryos. Thus long-range propagation of ‘localized’ 
effects may represent an important and general mechanism in 
dynamically guiding the emergence of size-dependent features of 
patterning systems. 

Our current simulation studies were not designed to system- 
atically and quantitatively fit all parameters to experimental data, 
an analysis that will require future efforts and data. In addition, 
all of the parameters used in our current simulation studies had 
been derived from fitting to training data sets that did not have 
the enhanced spread in embryo size. Naturally our simulated 
results are not expected to match perfectly with experimental 
observations for all the individual aspects even on a qualitative 
basis. In simulated large embryos, both the Hb2 and Krl 
boundaries shift posteriorly (Fig. 6a-h), whereas either an 
anterior shift or no significant shift is observed at early and later 
times of ncl4, respectively, in the experimental data (Fig. 2k]; 
Supplementary Fig. 2). The posterior Hb domain is expressed at 
lower levels in simulations compared with experimental data. 
Both of these disagreements with data occur at the edges of the 
modelled region, and might stem from regulators omitted in the 
model—the head gap genes and huckebein in the anterior and 
posterior, respectively. In addition, both empirical°°°* and 
recent modelling studies?” suggest a complex regulatory 
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Figure 6 | Simulation of Bcd scaling in the gene circuit model of gap genes. (a-h) Model output for Hb, Kr, Gt and Kni during cleavage cycle 14, at 3.125 
(early T1), 9.375 (late T1), 15.625 (early T2), 21.875 (late T2), 28.125 (T3), 34.375 (T4), 40.625 (T6) and 46.875 (T8) minutes after the thirteenth 
nuclear division. Simulations using controll? or scaled Bcd are shown as solid or dashed lines, respectively. The scaling of Bcd in large embryos was 
based on data from Cheung et al.4© (see Methods for details). The control Bcd profile was the median profile from wild-type embryos used by 
Manu et al.'? to fit the model. The invariance of the Bcd length scale in absolute units was represented by reducing the length scale by 20%—from 
0.1651 EL to 0.1321 EL—in the relative position units utilized in the model. Furthermore, a stronger Bcd source, caused by the increased amount of 
bcd mRNA in large embryos“®, was modelled by multiplying the source strength (A), and hence Bcd concentration, by 1.6. These values of Bcd amplitude 
and length scale imply a critical position at € = 0.311 (see Supplementary Fig. 9 for simulation experiments of adjusting the critical position). Control and 
scaled simulations are indicated in the legend above a. (i) Time series of boundary positions in the model output; see above for time points. Anterior and 
posterior boundaries are represented by left- and right-pointing triangles, respectively. The colour key is shown in b. Control and scaled Bcd simulations 
are shown as open and filled triangles, respectively. (j) Time series of the expression level (in a.u.) at the peak of the posterior Hb, central Kr, posterior 
Gt and abdominal Kni domains. Control and scaled simulations are shown as solid and dashed lines, respectively. Kr2, Kni2, Kni3 and Gt5 boundaries 
shift to the anterior. The reduction in Kr expression at the domain peak is relatively constant in time, whereas the difference between control and scaled 
Bcd simulations for Gt and Kni increases in time. 
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relationship between Hb and Kr in forming the Hb2 and Krl 
boundaries. 

The value of our current simulation studies lies in providing 
a fresh conceptual framework that bridges two disparate 
experimental observations—the long-range effects of low Kr 
expression and the under-scaling of Bcd in the posterior— 
through the well-documented and modelled regulatory properties 
of the gap gene network to explain the emergence of size- 
dependent features in the posterior. Our simulation results show 
that a reduced Kr level can on its own exert long-range and 
dynamic impacts on gap expression properties (in terms of both 
expression level and boundary position) in the posterior. The 
simulated long-range impact of Kr level reduction is reminiscent 
of the long-range effects on eve patterning observed in Kr and 
Kr/ +, but not kni, mutant embryos®!. A gene circuit model for 
Kr mutant embryos”? inferred that the anterior shift of the Gt 
boundaries was caused by a lack of Kr repression; our simulations 
suggest that a ~1.6 dosage of the Kr gene has similar effects, 
albeit of modest size. It is notable that the gene circuit model 
exhibits the Kr-mediated effects of size-dependent Bcd alteration 
without any additional training. Additional training with 
expression data from the size-selected lines would probably 
change the genetic interaction matrix, a zygotic property. Our 
simulations therefore suggest that the patterning differences 
originate maternally and do not require any divergence of the 
zygotic gap genes between the two lines. A related final point 
concerns the ability of the gene circuit model to both canalize 
high Bcd variation?)°° and exhibit posterior under-scaling as 
reported here. These apparently contradictory behaviours may be 
reconciled by noting that the perturbations are qualitatively 
different. The simulations in Manu et al. correspond to random 
Bcd variation in nominally isogenic embryos, whereas the 
perturbation modelled here is a systematic difference stemming 
from a large change in embryo length between two genetically 
distinct lines?!. That the same gene circuit underlies both types of 
the simulated—and empirically observed—behaviours highlights 
the complexity and plasticity of gene regulatory network 
dynamics under different perturbations. 

An important suggestion of our current study is that the Kr 
mRNA level divergence in relation to size can stem from the 
scaling properties of the Bcd gradient at the location of the Kr 
expression domain. Thus the AP patterning behaviour in the 
posterior part of the embryo with respect to scaling could be 
viewed/explained as being controlled primarily by the dynamic 
operation of the gap gene network, where the system’s initial state 
is determined by the Bcd input and mediated by Kr. The role of 
Bcd as a positive input for Kr expression is well documented both 
experimentally and in silico!??>*8-°°, and our current simulation 
results further support this relationship. Importantly, our study 
suggests that this relationship has a role in mediating the long- 
range and dynamic effects on gap gene expression to guide the 
dynamic emergence of size-dependent features in the posterior 
trunk region of the embryo. This supports the view that the entire 
AP axis is patterned as a single unified dynamic system in which 
the size-driven changes in maternal positional information are 
interpreted locally in individual nuclei, but in a globally coordinated 
manner. It remains to be elucidated whether our suggested Bcd-Kr 
relationship is the sole mediator for size-dependent patterning 
behaviour in the posterior at a quantitative level. 

Our experimental design takes advantage of two inbred lines 
that were selected from wild populations for laying large and 
small eggs*’. The fact that the spatial registry of gap expression 
boundaries and their temporal movements in these lines exhibit 
an overall similarity to one another and to other well- 
characterized wild-type strains suggests that neither of these 
lines is genetically ‘defective’ with respect to AP patterning. To 


further evaluate whether our observed Kr expression level 
divergence is indeed a property that is associated with embryo 
size (as opposed to other genetic differences between the inbred 
lines used in this study), we analysed Kr expression level in 
relation to embryo length within individual lines separately. Our 
analysis reveals a significant (inverse) correlation in embryos 
from the large-egg line, though not in embryos from the small- 
egg line (Supplementary Fig. 10a,b). These results document that 
Kr level can diverge in embryos where size differences are not 
caused by genetic differences. These results also further support 
our suggestion that Bcd gradient properties associated with 
embryo size, that is, its scaling properties, can impact Kr 
expression level. We note that, for extracting scaling-specific 
features, within-line analyses are inherently less sensitive than 
between-line analyses due to the limited size variations within a 
line (Supplementary Fig. 10). We attribute our inability to detect a 
significant correlation between Kr level and embryo length in the 
small-egg line to this reduced experimental sensitivity. This also 
further underscores the value of the experimental strategy used in 
our current study. 

Developmental robustness is an important concept that can 
guide how we view and interpret the relevant events that take 
place on different time scales. It has been proposed that within- 
species scaling properties of the Bcd gradient can be traced 
fundamentally to a dynamic relationship between bcd gene copy 
number expansion and nurse cell size expansion during 
oogenesis®. Under the framework of a recently proposed model 
of Tissue Expansion-Modulated Maternal Morphogen Scaling 
(TEM?S)°, this dynamic relationship also imposes constraints on 
Bcd gradient properties in relation to embryo size in different 
dipteran species during evolution. Relative to D. melanogaster 
embryos, the larger-sized D. yakuba embryos appear to exhibit 
both weakened gap expression levels (for Kr and posterior kni 
and gt domains) and under-scaled characteristics in the 
posterior’. This suggests a possibility that scaling-specific 
features and their mechanistic underpinnings within a species 
may also be applicable to evaluations of AP patterning properties 
across species. The early Drosophila embryo undergoes rapid 
nuclear divisions and morphological changes?*+416263. How cells 
(nuclei) in this dynamic environment decode positional 
information remains an open question??3!32346465_ The 
observed wave-like initiation and recovery of under-scaling in 
the posterior trunk region suggest that, to fully understand 
developmental robustness with respect to scaling, knowledge of 
the process of achieving the final outcome (that is, phenotype) is 
as important as that of the outcome itself. In addition, the 
temporal constraints on molecular processes imposed by the 
rapid morphological progression of the embryo may also be 
relevant to events and changes that take place on the evolutionary 
time scale. 


Methods 

Embryo collection and RNA FISH. The two inbred lines (#2.46.4 and #9.17.1) 
used in this study were gifts of Drs Cecelia Miles and Martin Kreitman and they 
had been selected from wild populations for laying large or small eggs*?. Flies were 
reared at 25°C and 0-4-h embryos were collected from 5-10-day-old females. For 
convenience, we refer to the embryos collected from the two inbred lines as large 
and small embryos, respectively. For mRNA FISH, digoxigenin (dig)-labelled RNA 
probes were synthesized in vitro using PCR products or plasmid (for hb) as 
templates as before®®. Gene-specific PCR primers are listed in Supplementary 
Table 2. The primers for kni (forward: 5’-CGGAATTCAGTCCTTCTTTGGCCGC 
TCTTAC-3’; reverse: 5/-TAATACGACTCACTATAGGGTCTAGATTAGACA 
CACACGAATATTC-3’) and the plasmid for synthesizing hb probes (CZ4001) 
were as described®. In our mRNA FISH analysis (refs 33,66), the hybridization 
step was carried out for 40h at 60°C. After hybridization with dig-labelled RNA 
probes, embryos were incubated overnight at 4°C with sheep anti-dig primary 
antibody (Roche, 1:400), followed by a 1-hour incubation at room temperature 
with goat anti-sheep secondary antibody conjugated with Alexa Fluor 594 (Life 
technologies, 1:400). DNA counterstain with 4,6-diamidino-2-phenylindole 
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(DAPI) (Sigma, 1:1,250) was performed for 10 min at room temperature before 
mounting using the Prolong Gold medium (Invitrogen). Coverslip ‘bridges’ were 
used to minimize physical distortions and the mounted embryos on slides were let 
to equilibrate for at least 10h at 4°C before imaging. 


Image acquisition. To allow comparisons of the expression profiles of a gene 
between embryos from the two inbred lines, we ensured that embryo collection, 
mRNA FISH and imaging were all performed on a side-by-side basis. During 
imaging, selection criteria were used to ensure that an embryo to be imaged: 

(1) was at the stage of interest (nc13 or ncl4 before gastrulation); (2) was judged to 
be laterally oriented with the anterior and posterior poles aligned horizontally 
(that is, without obvious tilt along the AP axis); and (3) did not exhibit gross 
deformations due to experimental manipulations or morphological defects. 
Image acquisition was performed on Zeiss Imager Z1 ApoTome microscope in 
conjunction with the AxioVision 4.8 software. All images were captured on the 
mid-sagittal plane of the embryo under linear settings without any normalizations 
or other nonlinear adjustments*!. To avoid pixel intensity saturation, we first 
selected a stained embryo that was judged to have the highest fluorescence signals 
(according to software determination) for calibrations. In particular, we used this 
embryo to set an exposure time such that the highest intensity pixels remained 
unsaturated, and this chosen exposure time, together with all the other microscopic 
and software settings, was fixed throughout an entire imaging cycle for all images. 
Therefore, absolute intensity levels for each gene could be compared between 
embryos from the two inbred lines, but not between genes irrespective of line. For 
each embryo, two images were taken, the first of which captured the fluorescent 
signals ( x 10 objective, capturing the entire embryo at its mid-sagittal section), 
followed by that for the DIC channel ( x 20 objective, focused on the dorsal side for 
quantifying the progression of membrane invagination). All images (in TIF format) 
had 1382 x 1040 resolution and 8-bit pixel depth. 


Time class classification. The nuclear cycle (nc13 or nc14) of an embryo was first 
determined by evaluating the number of nuclei on the dorsal sidef’. For embryos at 
ncl4 before gastrulation (~ 60 min into the interphase), we divided them into nine 
time classes based on nuclear morphology and membrane invagination. Nuclear 
morphology provides a sensitive measure of time for embryos at early times of 
ncl4 (refs 63,67), whereas membrane invagination is a reliable measure of later 
time?®296869, In our analysis, the nuclear length of an embryo was obtained 

by averaging measurements of five nuclei on the dorsal side captured by the 
DAPI channel. Membrane invagination ratio was calculated as a percentage of 
invagination depth relative to cortex length. We made measurements at three 
different locations on the dorsal side of an embryo and used the average for 
ranking embryos. Distance measurements for both nuclear length and membrane 
invagination depth were conducted using Image]. 

Embryos that had <15% membrane invagination ratio on the dorsal side, 
estimated not to have exceeded ~ 25 min into the interphase”, were grouped into 
the early time classes T1 and T2, whose division was based on the simple ranking of 
the measured nuclear lengths. Embryos that had > 15% membrane invagination 
ratio were divided into the remainder of the time classes (T3-T9) according to the 
simple ranking of the percentage of membrane imagination at arbitrary cut-offs 
(25, 35, 45, 55, 65 and 75%). The chosen method led to both a fine temporal 
resolution and a sufficient number of embryos at individual time classes for 
each of the two inbred lines and for each of the genes analysed (see Supplementary 
Table 1). We estimated that the time classes of T3-T5 represent mid-ncl4 embryos 
(~25-40 min into ncl4) and those of T6-T9 represent late ncl4 embryos 
(~40-60 min into nc14)®*. The nine time classes have an approximate time of: 
0-13.5, 13.5-25, 25-31, 31-35, 35-38, 38-42, 42-45, 45-48 and 48-60 min after 
the thirteenth nuclear division, respectively (see Fig. 6 for further details). 


Image processing to obtain gene expression profiles. Image processing was 
performed using Matlab (MathWorks) as follows*!. Briefly, an embryo was first 
oriented along its AP axis with the anterior end to the left and dorsal side up using 
the software with manual supervision when applicable. In our analysis, the length 
of an embryo was first recorded and divided into 50 equal-sized bins along the AP 
axis; this allows the extraction of a bin-location-based intensity profile on the 
dorsal side of the embryo such that the extracted profile has effectively been 
‘projected’ onto the AP axis and expressed as a function of relative AP position 
irrespective of embryo size or shape. Specifically, the mRNA FISH intensities from 
the image of the embryo were extracted by moving a scanning circle (of 61 pixels) 
along the dorsal side of an embryo immediately below (basal to) the nuclear layer 
as detected by DAPI signals. Background signals on the DAPI channel were also 
used for defining the embryo bounds with manual assistance whenever necessary. 
Intensity values for each scanning location were expressed, in arbitrary units, as the 
mean of all pixels within the scanning circle and an embryo’s raw intensity profile 
smoothed by applying a Gaussian filter. Fluorescent signals at non-expression 
regions were then used in an embryo-specific background subtraction across all 
positions, without any further adjustments unless otherwise specified. 


Boundary position and movement and domain width. To determine the 
boundary position of a gap expression domain, the mean of three highest intensity 


values of a domain was set as 1 for normalizing intensity values at all positions for 
this domain. Boundary position was defined as the relative AP position (€ = x/L) at 
which the normalized intensity value is 0.5 (see Fig. 4a). For the posterior tll 
domain, we avoided the most posterior three bins for normalization and boundary 
calculation to prevent capturing erroneous intensities from pole cells. We also 
avoided boundary calculations whenever prohibited because of low mRNA 
abundance that rendered such measurements unreliable. In particular, the 
boundaries of gt] and gt2 become detectable only after T4, whereas the mRNA level 
of the anterior hb expression domain decreases so rapidly during early ncl4 that 
the anterior boundary (hb1) no longer exists after T4 (Supplementary Fig. 1). Our 
analysis also did not include mRNA expression domains that never reached an 
appreciable level suitable for reliable boundary measurements (such as the 
posterior Kr domain). Thus the data shown in Fig. 1c,d and Supplementary Fig. 2 
do not have boundaries for certain genes at certain times. 

The width of an expression domain was defined (see Fig. 4a) as the distance 
between its anterior and posterior boundaries (both expressed as relative AP 
positions). Specifically, for each time class, we first measured the domain width of 
an expression domain in individual embryos to obtain a mean (+ s.d.) of the large 
or small embryos as a group, (Wiarge) and (Wemai), respectively. We then 
calculated the width difference at each time class as AW = (Wiarge) — (Wsman)- The 
s.d. of AW was estimated by \/o7/m + 03/n2, where n, and n, are sample 
numbers for the two groups of embryos, and a, and o are their respective s.d. The 
gt ant2 domain is composed of two sub-domains upon its split at later nc14. For the 
posterior tll domain, the domain width is defined as the distance between the tll3 
boundary and the posterior end. When calculating boundary movements as a 
function of time, s.d. of the difference (distance) of the means was also estimated by 


Voi /m +03/m. 


Scaling coefficient and threshold-crossing analysis. To analyse the scaling 
property of an expression boundary, we pooled the large and small embryos at a 
given time class to generate a scatter plot of boundary position (č) against nor- 
malized embryo length L/(L), where (L) is the mean of the two inbred lines’ mean 
embryo lengths in this study. The slope of the fitted linear regression in such a 
scatter plot is defined as S (see Supplementary Fig. 3). Threshold-crossing analysis 
was performed using the interpolated mean (unnormalized) expression profiles of a 
group of embryos to obtain positions at arbitrary thresholds. Data analyses and 
figure generation were performed using Matlab (MathWorks) or Microsoft Excel. 
Unless otherwise specified, all values shown are mean + s.d. and all P values are 
from two-tailed Student’s t-test. 


Simulations. The simulations were performed using the gene circuit model of 
Manu et al.!947, The model simulates the gene regulation of Hb, Kr, Gt and Kni in 
nuclei lying between 0.35 and 0.92 EL on the AP axis of the embryo. The 
simulation begins after the twelfth nuclear division and lasts until gastrulation, a 
duration of ~71 min (ref. 62). The state variables are the concentrations of each 
gap protein in each nucleus, v?, where a is the gap gene and i is the nucleus. The 
concentrations evolve according to the coupled differential equations 


dv? i 
m =R°g (= Tye mpi 5 E” wh (t) — r) a) 
b b 


+D. (vf, Hii — 2v/) — & v7, 


where the first term represents synthesis, the second represents the diffusion of gap 
proteins, and the third represents first-order degradation. R°, D and 6° are the 
maximum synthesis rate, diffusion coefficient, and degradation coefficient, 
respectively. In the synthesis term g(u) = 1] matt i is a sigmoidal 
thresholding function. T is the genetic interconnectivity matrix between the gap 
genes, in which positive and negative elements signify activation and repression, 
respectively. v>“¢ is the concentration of Bcd and m* is the strength of regulation by 
Bcd. The model also includes the time-dependent effects of upstream regulators 
Cad and TIl, where w?(t) is the concentration and E?” is the strength of regulation. 
h’ is a threshold term. The initial conditions for the gap genes were from nc12 data 
in the FlyEx database (http://urchin.spbcas.ru/flyex/)°>”°. Kr, Gt and Kni protein is 
first detected in nc13 (ref. 29), and hence their initial concentrations were zero. Cad 
and TII data were from?!. Except as noted below, we utilized the model reported in 
Manu et al.!? (model 007) without modification, and refer the readers to that paper 
for the details of numerical simulation and parameter inference. The parameter 
values are listed in Supplementary Table 3. 

We designate the unmodified model, which utilizes a Bed profile vB4(x) = A 
exp( — x/2) with A = 351 (in a.u.) and À = 0.1651 EL, as the control simulation. For 
the simulation of reduced synthesis rates (Supplementary Fig. 8), one of R, REY, 
R@, RK" or all together were reduced by 20%. For the simulations of Bcd scaling in 
larger embryos (Fig. 6 and Supplementary Fig. 9), we considered two effects. First, 
the length scale of Bcd does not change with embryo length*®, a fact we represent 
by reducing the relative-unit length scale of Bcd by 20%—the approximate ratio of 
embryo lengths of large and small embryos—to 0.1321 EL in the model. Second, 
larger embryos were observed to have more bcd mRNA and higher protein 
expression at the anterior*®. This was represented in the model by increasing the 
value of A. The simultaneous reduction of relative length scale and increase in A 
causes the scaled Bcd to intersect with the control Bcd at a critical position“. We 
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investigated the dependence of simulated gap gene expression patterns on the 
location of the critical position by varying A (Supplementary Fig. 9). Larger values 
of A shift the critical position toward the posterior and vice versa. 

We measured the boundary positions and peak expression level in domains by 
fitting a smoothed cubic spline to the model output using the CSAPS function of 
MATLAB (MathWorks Inc.). Boundary positions were calculated as the position 
where the simulated expression level is at the mean of the consecutive extrema. 
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